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Stable vortex iV-omers are constructed in coherently coupled iV-component Bose-Einstein conden- 
sates. We classify all possible iV-omers in terms of the mathematical graph theory and numerically 
construct all graphs for N = 2,3,4,5. We also find that iV-omers are well described as CP N ~ X 
skyrmions when inter-component and intra-component couplings are U(N) symmetric, and we eval- 
uate their size dependence on the Rabi coupling. 
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Introduction — Exotic vortices in multi-component 
condensates have recently attracted considerable atten- 
tion in the area of condensed matter physics; in this 
light, the topics of interest include exotic superconduc- 
tors, superfiuid 3 He, Bose-Einstein condensates (BECs), 
exciton-polariton condensates, and nonlinear optics. In 
particular, recent advances in realizing BECs in ultra- 
cold atomic gases have opened up new possibilities in 
quantum physics [H [2] ; of particular interest is the tun- 
ing of the s-wave scattering wavelength via a Feshbach 
resonance [3H5]. Multi-component BECs can be real- 
ized when more than one hypcrfine spin state is si- 
multaneously populated E] or when more than one 
species of atoms are mixed [H [7J [8] . The recent exper- 
imental achievement of a condensate of ytterbium of- 
fers condensations up to five components [5J. Vortices 
in multi-component BECs have been realized experimen- 
tally [TU1 |TT] , and the structures of these vortices are con- 
siderably richer than those formed of single components 
P2H2S] ; one fascinating feature in such cases is that vor- 
tices are fractionally quantized. The study of BECs can 
provide an ideal opportunity to examine the dynamics of 
exotic vortices because these vortices can theoretically be 
explained via a quantitative description using the mean 
field theory by the Gross-Pitaevski (GP) equation and 
further, the vortices can be experimentally controlled to 
a large degree. 

Vortices or fluxes are fractionally quantized in multi- 
gap superconductors [26-29 as in the case of multi- 
component BECs. One added feature of multi-gap su- 
perconductors is the presence of Josephson coupling be- 
tween the gaps. Due to the presence of this coupling there 
appears a sine-Gordon soliton for the phase difference be- 
tween the two gaps .26[ [3D]. When a set of two vortices 
in two different gaps is separated, a sine-Gordon domain 
wall connects the vortices, thereby resulting in a two- 
vortex molecule, i.e., a dimer |29| . Multi-components 
are also present in thin-layer Josephson junctions and 
high-T c superconductors with multi-layers. Recently dis- 
covered iron-based superconductors may have three gaps, 



in which case there may be three- vortex molecule i.e., a 
trimer [3T]. However, the vortex molecules have not been 
observed directly thus far in superconductors except in- 
directly [32], because the repulsion between vortices is 
exponentially suppressed by the Meissner effect so that 
it cannot be balanced with the domain wall tension. 

In the case of multiple hyperfine spin states of BECs, 
internal coherent (Rabi) couplings between multiple com- 
ponents can be introduced by Rabi oscillations, similar to 
Josephson couplings in multi-gap superconductors. As in 
superconductors, a sine-Gordon domain wall of the phase 
difference of two components appears |33| . In this light, 
the advantages of BECs are that Rabi coupling can be 
tuned experimentally and that a vortex dimer can exist 
stably because repulsion between two fractional vortices 
[20] is sufficiently strong in the absence of the Meissner 
effect, to be balanced with the tension of a sine-Gordon 
domain wall connecting them |17j . Recently, we observed 
stable-vortex trimers in three-component BECs, and we 
showed that the shape and the size of the molecule can be 
controlled by varying he strength of the Rabi couplings 

ED 

In this Letter, we construct stable-vortex ./V-omers, 
that is, molecules made of N fractional vortices in TV- 
component BECs, and we find that 7V-omers with N > 3 
exhibit several novel properties that dimers do not pos- 
sess, i.e., the existence of chirality pairs and metastable 
states such as twist, holding, and capture, which are 
properties exhibited by chemical molecules. Each con- 
densate wave function has a nontrivial winding number 
around one of the N fractional vortices. When the Rabi 
couplings are turned on, the fractional vortices are glued 
by domain walls to form a vortex A^-omer. We make 
use of the mathematical graph theory to classify the vor- 
tex A^-omers whose number exponentially increases as N 
increases. We numerically construct all possible graphs 
for A~ = 3,4, 5 by imaginary time propagation with the 
phase winding and a constant density fixed at the bound- 
aries, thereby our results imply that the vortex Af-omers 
are formed in rotating BECs. We find that two geomet- 
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rically different graphs topologically equivalent to each 
other can both be stable. Furthermore, we can con- 
trol the shape and size of the molecules by varying the 
strength of the Rabi couplings. To demonstrate this, we 
construct a vortex constellation, i. e., the orion. We also 
show that vortex TV-omers are well described by CP^ -1 
skyrmions when all possible Rabi couplings are set to the 
same values, where CP^" 1 ~ SU(N)/[SU{N-l)xU{l)} 
denotes the complex projective space. 

TV-component BECs — We consider coherently cou- 
pled TV-component BECs of atoms with mass m, de- 
scribed by the condensate wave functions 'I'; (i — 
1, 2, • • • , TV) with the GP energy functional 

i,j=l J \ m 

-/*i|*i| 2 <5 ii - w «*i*A (!) 

where atom-atom interactions are characterized by the 
coupling constants gij = 4ttH aij/m with s-wave scatter- 
ing lengths dij, fii denotes the chemical potential, and 
LJij = ujji [uja = 0) denotes the internal coherent cou- 
pling due to Rabi oscillations between the i-th and j'-th 
components. We mainly consider the case where gu = g, 
and fa = jU, and g^ = g (i ^ j) in this Letter. It is 
straightforward to consider the general case but all the 
results below are essentially unchanged from our simplest 
choice. 

The symmetry of the Hamiltonian depends on the cou- 
pling constants g,g, and LJij. The maximal symmetry 
U(N) is achieved when g = g with = 0. When g ^ g 
and LJij = 0, this symmetry reduces to U(1) N acting on 
a phase of each component. When g ^ g and Wy ^ 0, 
it further reduces to the gauge symmetry 1/(1) which ro- 
tates the phases of all the components uniformly. 

When all the internal coherent couplings are equal i.e., 
uj.ij = uj, the condensations of the ground state are — 

v = J g^jvlj j^; w ith i = 1, 2, • • • , TV. These amplitudes 
are modified in the case LJij ^ utyj', and one should solve 
the variational problem SE/d^i = to determine Uj = 
l^/jj numerically. 

As long as the internal coherent couplings |cay| are 
maintained sufficiently small with respect to the other 
couplings (we choose g ~ g = O(10 3 ), p = O(10 2 ), and 
u = O(10~ 2 )), the symmetries U(N) for g = g and 
U(1) N for g 7^ g are nearly intact. These symmetries 
are spontaneously broken in the ground state. Conse- 
quently, vortices that are quite different for g = g and 
g > g can appear. For g = g, an axisymmetric giant 
vortex appears. As we will explain below, this can be 
interpreted as a CP N_1 skyrmion. On the other hand, 
for g > g, there appear TV fractional vortices associated 
with the broken U(1) N symmetry, which are connected 
by domain walls, thereby resulting in molecules with N 



vortices. 

The internal coherent couplings — u>ij^*^ j = 
—2viVjU>ij cos(6i — 6j), with 6i — arg'I'i, glue the frac- 
tional vortices by creating domain walls among them. For 
LJij > 0, all the phases 9i and 6j coincide, i.e., 6{ = 0j in 
the ground state. 

First, let us consider the case that all Wy's are zero. 
A vortex winding around only one component, ^ = 
Vip(r)e , in polar coordinates r and 9, with the other 
components almost constant, is fractionally quantized as 

/v 2 h 
dr-v s = ^-, (2) 

with the superffuid velocity v s , and it is stable. How- 
ever, once we turn on the internal coherent couplings 
lj^, no matter how tiny they are, the vortex is in- 
evitably attached by semi-infmitely long domain walls 
extending to the boundary because of the Rabi energy. 
To examine this phenomenon, we approximate the con- 
densate wave functions by = Vie l0i to obtain the 

reduced energy functional £ ph asc = J2iLi l^(^) 2 ~ 
^■Y^i>j UJ ij v 'i v j cos(#i — 9j). This reduced energy func- 
tional shows that TV — 1 domain walls are attached to 
the i-th vortex in the direction Qi — 8j. Therefore, this 
vortex cannot remain stable in an isolated state, and it 
must be connected to the boundary or to vortices wind- 
ing around the other components forming the molecules. 

As described below, we find that the mathematical 
graph theory is useful to classify the vortex molecules, 
where we identify vertices as vortices, and the edge con- 
necting i-th and j-th vertices as the Rabi coupling u>ij if 
it is nonzero. 

<CP N ~ 1 skyrmions — Here, we consider the symmetric 
case with g = g and with all possible internal coherent 
couplings being the same (aty = lo for all which cor- 
respond to the complete graphs. We numerically solve 
the GP equation by using an imaginary time propaga- 
tion for TV = 3,4,5, as shown in Fig. [I] We find that 
energy densities are always axisymmetric so that we may 
call them giant vortices. However, each component has 
a different configuration; the vortex centers of each com- 
ponent are expressed as dots, whose positions are Z^v 
symmetric. We also varied w and found that the size of 
the configuration decreases as lo increases, as shown in 
Fig. [2] for TV = 3,4,5. With our numerical survey over 
a wide range of parameter space and various initial con- 
figurations, we conclude that the existence of vortices is 
robust. 

We further confirm the stability and existence of the 
solutions by using a variational method. In the parame- 
ter region of g 3> |cd|, we note from Eq. ([!]) that the total 
density is approximately constant, i.e., ^\ \^i\ 2 ~ p/g- 
Thus, we are naturally led to define the reduced wave 
functions {fa} by *j = <T> 1; J^iLi \<t>i\ 2 = !• Con- 

sequently, the GP model expressed by Eq. (TH) reduces to 
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FIG. 1: Energy-density distribution of axisymmetric 
molecules with m — h — 1, g — g = 10 3 , fi = 10 2 , and 
cj = 0.05. The size of the simulation box is (2L) 2 = 20 2 with 
the lattice space set to be less than 0.1. Zeros of are indi- 
cated by the black dots. The numerical range for each of the 
figures is set as x G [—4,4] and y £ [—4,4]. 
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FIG. 2: Sizes of CP N ~ 1 skyrmions for N = 3,4,5. The dots 
indicate skyrmion radii obtained from the numerical solutions. 
The radii obtained by the variational method are connected 
by broken lines. 



the S nonlinear si^ma model. Upon separating the 
total phase 9 = ^fYli=i^i> the rest of the degrees of 
freedom are described by the CP^" 1 model, 



-CP" 



JL 

Ng 



2m 



,(3) 



with the identification fa ~ e %a fa. Motivated by exact 
solutions of CP N ~ X skyrmions and the In symmetric 
zeros in Fig. [I] we propose a variational ansatz 
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with variational parameters a and A, as a generalization 
of the case of N = 2 (CP 1 ) [H [IS]. Substituting the 
above values into Eq. ([3|, we find 



NgE 

2"K[L 



dr r 



h 2 r 2 + X 2 e- ar ' 2 + (ar 2 



2m 



r 2 + A 2 e 



A 2 £ 



-(N-l) 



(6) 



where L denotes the system size. We next find the values 
of A and a for minimizing the energy and read the size 
from r = Ae~2 r . In our study, we observed complete 
agreement between the radii obtained by the numerical 
solutions and the variational method, as shown in Fig. [2] 

Vortex molecules — Next, let us consider the miscible 
case g > g to observe non-axisymmetric configurations, 
namely, the vortex molecules. 

Dimers: Vortex dimers are molecules of two vortices 
appearing in two-component BECs [17] . The dimer is the 
simplest molecule that does not have any degeneracies or 
substructures. 

Trimers: Vortex trimers are molecules of three vor- 
tices appearing in three-component BECs. These exhibit 
two new properties that the dimers do not have. The 
first reported in Ref. [23] is that the shape of the trian- 
gle can be changed by tuning the parameters in Eq. ([I]). 
The second is the chirality of the triangle. This can be 
easily seen when the coupling constants are generic, so 
that the vortices have different sizes, as seen in Fig. [3] 
In two-component BECs, chirality does not make sense 




FIG. 3: Chirality: two trimers in a three-component BEC 
(flu 7^ ff22 7^ 333 7^ <?ii)- The plotted energy densities are 
shown in the figures. In each figure, there are three different 
vortices forming a molecule. The left panel is a mirror image 
of the right panel. 



because left and right can be exchanged by tt rotation. 
On the other hand, as shown in Fig. [3j the left and right 
configurations cannot be transformed to each other by a 
rotation but can be transformed by a mirror. The trimers 
are energetically completely degenerate. 

Tetramers: In multi-component BECs with more 
than four condensate wave functions, the number of pos- 
sible molecules increases exponentially. Therefore, in or- 
der to classify such vortex molecules, we make use of 
the mathematical graph theory, in which vortices are ex- 
pressed by vertices and the Rabi couplings are expressed 
by edges. Graphs isomorphic to each other are not distin- 
guished in the graph theory. The number of independent 
connected graphs relevant for four-component BECs is 
six. In the graph theory, a graph is characterized by the 
sequence of the number of edges connected to each ver- 
tex. For instance, (1, 1, 2, 2) implies that two vertices are 
connected by two edges, and other two vertices are con- 
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FIG. 4: (a)-(h) All possible tetramers in a four-component BEC and (i) the orion in a seven-component BEC. We set 
g — 900 while the other parameters are the same as those for Fig. [T] For the tetramers, we set = 0.05 and 
uiij — for the connected and disconnected pairs, respectively. For the orion, we set (wi2, W23, ^34, W45, i^ie, <^56, ^37, ^67) = 
(0.02, 0.001, 0.02, 0.02, 0.05, 0.02, 0.2, 0.2) and the rests to be zero. The dots and dotted lines indicate the vortex center of each 
component and the Rabi couplings between components, respectively. The color maps represent the energy density (the left 
for tetramers and the right for the orion). The number at the upper-right corner in each subfigure represents the total energy. 
As regards more precise information on the tetramers, the phase and amplitude of each component are plotted in Figs. [5] and 
[5J and phase differences between two among four components are plotted in Fig. [7] in the appendix. In Fig. [7[ the positions of 
domain walls are visible. 



nected by one edge. We obtained numerical solutions for 
a four-component BEC with the same values of Rabi cou- 
plings (Fig. [4]) , exhausting all possible graphs with four 
vortices as vertices including six connected graphs seen in 
Fig. [4] (a)-(f) and two disconnected graphs in Fig. [3] (g) 
and (h) . We calculated the GP energy ([I]) subtracted by 
the ground state energy for each configuration, as shown 
at the top-right of each panel in Fig. |4j where we chose 
r = 6 for the spatial integration in Eq. ([T]). 

Our simulation indicate the occurrence of several new 
phenomena that do not exist in dimers or trimers: 1) 
twist, 2) holding, and 3) capture. We have found that 
most molecules are accompanied by "twisted" molecules, 
as in Fig. [4] (a)-(e). A pair of untwisted and twisted 
molecules corresponds to identical graphs isomorphic to 
each other, while both pairs are energetically stable, cor- 
responding to absolute and local minima; the upper con- 
figuration has lesser energy than the lower configuration 



in each figure in Fig.|4](a)-(e), and corresponds to the ab- 
solute minimum. There exists an energy barrier between 
twisted and untwisted molecules, and therefore, once a 
twisted vortex molecule is formed, a certain amount of 
finite energy is required to "untwist" it. As an example, 
a rod molecule denoted by the graph (1,1,2,2) in (a) 
has an energy value of 1.17 when it is straight, while its 
energy value is 1.21 when it is twisted. 

The second and third new phenomena can be observed 
in the disconnected graphs (g) and (h). The holding phe- 
nomenon can be observed in Fig. [4] (g) . This is seen as 
the union of the two graphs with (1,1). If the two domain 
walls connecting vortices are not crossed, this molecule 
breaks up into two dimers because of repulsion between 
the vortices. The third phenomenon observed is the ab- 
sorption of a molecule inside a bigger molecule, as in 
Fig.|(h). 

In the mathematical graph theory, an energy compo- 
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nent called the graph energy is defined motivated by the 
7r-electron energy of molecules in chemistry |35j . In such 
a case, for a given graph, configurations with the mini- 
mum graph energy are unique. Our GP energy ([!]) pro- 
vides a new definition of graph energy in the graph the- 
ory, which makes vortex molecules more interesting even 
mathematically. 

A-omers: We can engineer as many multiple vortex 
molecules as we want with the constituent vortices. As 
an example, we present vortex pentamers corresponding 
to all possible connected graphs with five vertices in a 
five-component BEC with the same Rabi couplings in 
Fig. [H] in the appendix. 

Thus far, we have concentrated on the case when 
ujij = lj 7^ and u)ij = for the connected and dis- 
connected pair of the vortices, respectively. We can con- 
trol the positions and shape of the molecule by varying 
tOij inhomogeneously. As an example, we present a vor- 
tex heptamer (seven vortices), designed as the orion in 
Fig.g](i). 

Conclusion — In conclusion, we have demonstrated 
the construction of stable-vortex A-omers in coherently 
coupled A-component BECs, and we have classified all 
possible A-omers in terms of graph theory. One graph 
may correspond to two stable configurations which are 
the local and global minima. We observed several new 
phenomena including chirality, twist, holding, and cap- 
ture. We also found that A-omers are well described by 
CP N ~ X skyrmions for the U(N) symmetric couplings. 
We used various initial conditions and always reached 
the same configurations (up to rotation), and therefore, 
we conclude that our results are consistent and verifi- 
able. The GP energy provides a new definition of graph 
energy in the graph theory. When the values of u)y 's are 
not equal, the domain walls have tensions different from 
each other and the corresponding edges have weights in 
the graph theory. When Uj's are not equal, the vortices 
have different masses and the corresponding vertices have 
weights. 

In experiments, an integer vortex is created as usual 
by gradually increasing the rotation, and subsequently, 
it is split into a vortex A-omer. Two-component BECs 
of different hyperfine states of the same atom have been 
already realized using the |l,-l) and |2, l) states [10] 
and the |2, l) and |2, 2) states [34] of 87 Rb, respectively. 
We believe that systems with three or more components 
can be realized by using a mixture of the above mentioned 
states of 87 Rb via an optical trap |36] , and our prediction 
is testable in laboratory experiments. 
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Appendix — Here, we provide more detailed informa- 
tion on our numerical solutions. Figures [5] and [6j show 
the plots of the amplitudes and the phases of condensate 
wave functions Vfj for the tetramers shown in Fig s. |4[ a)- 
(h) . We also show the relative phases of $j in Fig. [71 The 
figures clearly indicate that the condition 9i — 9j at the 
boundary is satisfied along with showing the positions of 
the sine- Gordon domain walls. 

Finally, we illustrate in Fig. [8] all the vortex pentamers 
in a five component BEC that are topologically incquiv- 
alent. 
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FIG. 5: The amplitudes and phases of V^i, ^2, ^3, and $4 for the tetramers shown in Fig. [4] The arrows indicate the phase 
vector (Re\&i, Im\&j). The colors indicate the amplitudes of the condensate wave functions. 




FIG. 6: See the caption of Fig. [5] 
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FIG. 7: The relative phases (arg*; -arg*j for {i,j} = {(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)}) for the tetramers shown in 
Fig. [4] The panels surrounded by the red squares correspond to the relative phases for the non-zero Rabi frequencies. Namely, 
they show the positions of the sine-Gordon domain walls. 
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FIG. 8: Representatives of all connected pentamers with the least energy in a five component BEC. All parameters are the 
same as those in Fig. [4j 



